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Abstract 


An approximate Riemann solver is developed for the governing equations of ideal 
m agnetohydro dynamics (MHD). The Riemann solver has an eight- wave structure, 
where seven of the waves are those used in previous work on upwind schemes for 
MHD, and the eighth wave is related to the divergence of the magnetic field. The 
structure of the eighth wave is not immediately obvious from the governing equations 
as they are usually written, but arises from a modification of the equations that is pre- 
sented in this paper. The addition of the eighth wave allows multi-dimensional MHD 
problems to be solved without the use of staggered grids or a projection scheme, one or 
the other of which was necessary in previous work on upwind schemes for MHD. A test 
problem made up of a shock tube with rotated initial conditions is solved to show that 
the two-dimensional code yields answers consistent with the one-dimensional methods 
developed previously. 


*This research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS 1-19480 while the author was in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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1 Introduction 


The governing equations of ideal magnetohydrodynamics (MHD) describe the physics of a 
conducting fluid in which the following assumptions hold: 
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where p, V, r and L are, respectively, characteristic density, speed, time and length scales 
for the problem, c is the speed of light, and e and a represent the dielectric constant and 
conductivity of the fluid. These equations, written in conservation-law form, are 
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where I is a 3 x 3 identity matrix, p is the density, u is the velocity, p is the pressure, B is 
the magnetic field, and E is the energy, defined as 
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( 3 ) 


Solutions of these equations can yield insight into a number of problems governed by fluid- 
dynamic and electromagnetic effects. 

Much of the past work in solving these equations has bden based on Rusanov and Lax- 
Wendroff techniques. Only recently have authors begun to work on upwind schemes for 
solving these equations. In particular, Brio and Wu [2], Zachary and Colella [?], and Dai 
and Woodward [7] have done some of the early development of Riemann-solver-based schemes 
for the MHD equations. Their work has been based not on the system of eight conservation 
laws as written in Equation 2, but instead on the closely related system that comes from 
assuming B x = constant and dropping the evolution equation for B x . This yields a 7 x 7 
system. The reason for their use of this modified system arises from the fact that one of the 
equations governing the magnetic field is 


V • B = 0 , 


( 4 ) 




Figure 1: Waves in the One-Dimensional MHD Riemann Problem 


which, in one dimension, becomes the constraint B x = constant. 

The eigenvalues and eigenvectors of this 7x7 system are well known (see, for example, 
the book by Jeffrey and Taniuti [3]); they correspond to: 

• one entropy wave traveling with speed u\ 


• two Alfven waves traveling with speed u ± c a where 



is the Alfven speed; 


• four magneto- acoustic waves, two “fast” and two “slow”, traveling with speed u ± cj 
and u ± c s respectively, where 

2 _ 1 ( 7P + B • B _L (lP + B • B^ 2 

j ~ 4 irJ ■ 

An (x,t) diagram of the wave interactions at a cell interface is shown in Figure I. 

Given these seven eigenvalues and corresponding right and left eigenvectors, it is pos- 
sible to develop a linear approximate Riemann solver ala Roe [2, ?], or a more nonlinear 
approximate Riemann solver [7]. Once some questions as to how to scale the left and right 
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eigenvectors of the system are answered (for a very nice description of the problems of scaling 
in the MHD eigensystem, and an elegant solution to these problems, see the paper by Roe 
and Balsara [5]), a robust solver for one-dimensional unsteady problems in MHD can be 
developed. 

Building a code capable of solving two- or three-dimensional problems from the one- 
dimensional Riemann solver building block is not, unfortunately, as straightforwaid as in 
the case of the Euler ecjuations. In the one-dimensional problem, no evolution equation 
is necessary for the component of B normal to a cell interface, because the condition of 
Equation 4 implies that B Ut = B UR . However, in a two-dimensional problem, this is no 
longer true. In two dimensions, the discrete constraint corresponding to Equation 4 is 

E £ n d.s = o, ( 5 ) 

faces 

and so a jump in B n is allowed across a face; it simply must be balanced by the jumps 
across the other faces of the cell. Thus, a separate procedure for updating this portion of 
the magnetic field must be implemented, and must be implemented in such a way as to 
satisfy the constraint implied in Equation 5. It should be noted that the VB constraint is a 
headache not just for upwind schemes for MHD, but for solution of MHD problems in more 
than one dimension by any method. Typically, one of three approaches is taken to satisfy 
this constraint: 

• a projection scheme, in which a Poisson equation must be solved to subtract off the 
portion of the magnetic field that leads to a non-zero divergence; 

• non-collocated variables (e.g. a staggered-grid approach), so that the constraint is met 
identically; 

• a vector-potential description of the magnetic field, so that the constraint is met iden- 
tically. 

A very different approach is taken in the work presented here. Instead of solving a seven- 
wave Riemann problem, with an added procedure to update the remaining B-field component 
which assures that Equation 4 is satisfied, an eight-wave Riemann solver, in which all of the 
magnetic field components are updated, is developed and tested. 

2 Derivation of the Eight- Wave Riemann Solver 

Given the primitive variables 

W = (p,u,v,w,B x ,B y ,B z ,p) , (6) 
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The Riemann solver would normally be based on the eigensystem of A p , but it is evident 
that this matrix is singular — the fifth row of the matrix is zero, leading to a zero eigenvalue. 
This zero eigenvalue is clearly non-physical — the eigenvalues should appear either singly as 
the x — component of the flow speed, u, or in pairs symmetric about u. It also does not bode 
well numerically — the mode corresponding to this eigenvalue will be undamped. 

The approach taken here is to look for a way in which to modify the governing equations 
so as to make A p non-singular. The criteria that should be met by the modified matrix A* 
are: 


• The eigenvalues and eigenvectors corresponding to the seven waves in the one-dimensional 
(B x — constant) Riemann solver remain unchanged; 

• The eigenvalue corresponding to the new eighth wave is equal to u (the only physical 
choice for a single eigenvalue); 

• The left and right eigenvectors corresponding to the new eight wave “make sense”; 

• In the case B x — constant, the eight- wave Riemann problem reduces to the seven- wave 
Riemann problem. 

With these criteria in mind, it becomes possible to find a modified version of A v , given 
some patience and some facility with Maple’s symbolic manipulation capabilities. The mod- 
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ified matrix that meets the above criteria is 


up 0 000 00 

0 u 0 0 0 ^ f J 

0 0 u 0 0 -&■ 0 0 

000 u 0 0 0 

00 0 0 u 0 00 

0 By -B x 0 0 u 0 0 

0 B z 0 -B x 0 0 « 0 

0 7 p 0 0 0 0 0 u 

The eigensystem of this matrix is composed of the following eight waves, with their corre- 
sponding eigenvalues A, left eigenvectors i and right eigenvectors f : 

One Entropy Wave 

A e = u 

l = (l,0, 0,0, 0,0,0,-^) 

f t = (1,0,0,0,0,0,0,0) T ; (10) 





Four Magneto-acoustic Waves 


A j, s —u± c/ tJ 
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One “Divergence” Wave 


A d = u 

Id = ( 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ) 

fd = (0,0, 0,0, 1 , 0, 0, 0) T . (13) 


It is important to note that the first seven waves yield the same eigenvectors and eigen- 
values as the seven-wave Riemann problem, with the additional information that none of 
them carries a change in B x (the fifth entry of each right eigenvector is zero), and none of 
the wave strengths is proportional to a jump in B x (the fifth entry of each left eigenvector is 
zero). The new eighth wave travels with the x — component of the flow speed (its eigenvalue 
is u), and it carries a jump in B x (the only non-zero entry in the left eigenvector is the entry 
corresponding to B x ), and affects only the x — component of the magnetic field (the only 
non-zero entry in the right eigenvector is the entry corresponding to B x ). 

It is clear that the eigensystem of this modified matrix has all of the desired properties. 
In the case B x — constant, the strength of the eighth wave is zero, and the model reverts to 
that of the seven-wave problem. The new wave simply gives a rational procedure for dealing 
with non-zero jumps in B x across the cell faces, which will in general occur when problems 
in two or three dimensions are being solved. The question remains, however, of what the 
modification of the matrix A p (and the corresponding changes to B v and C p ) has done to 
the system of conservation laws. 

This can be seen by collecting the source terms due to the modifications to A p , B p and 
C p and transforming to conserved variables. The new equation set, which has the eight-wave 
eigensystem described above, is 



P u 

puu + I (j) + — BnB 

uB — Bu 

(E + p+^)u-B(u-B) 



V B . 


(14) 


This is a noteworthy result: the source term that must be added to Equation 2 is proportional 
to V - B. At the partial differential equation level, only terms that are equal to zero have 
been added to the conservative form of the governing equations. So, while technically the 
equations are no longer in conservative form, the deviations from conservation will be very 
small. It is only by writing the equations in this slightly non-conservative form that the 
singularity related to V • B can be removed. It has been previously noted that solving the 
momentum equation in non-conservative form can remove instabilities related to non-zero 
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V B[l|; the current work hopefully reinforces this earlier result, and sheds further light on 
the mechanism for stabilizing the equations, as well as applying the idea in a novel way to 
develop a Riemann solver for multi-dimensional MHD. 

It is interesting to note another justification of this particular choice of source term. 
Rewriting Equation 2 slightly by expanding some of the terms, the following form of the 
equations 

% + V-(/m) = 0 
p + - B YB - BV B = 0 


(pu) + V ■ (puu) + V ( 


— + u YB + BY u - B Vu - uY B = 0 

dt 

OE B B\ / B B\ 

— +uVlE + p+ — 1 + IE + P+— JVu- 


B ■ V (u B) - (u • B) V - B 


(15) 


is obtained. The terms that are proportional to VB have been underlined; they are exactly 
the same as the source term defined above. Thus it can be seen that the addition of the 
source term in Equation 14 simply acts to remove the terms proportional to V • B that 
appear in Equation 2. 

Another interesting note is what the evolution equation for V • B is for the two forms of 
the governing equations. This may be seen by taking the divergence of the evolution equation 
for the magnetic field in Equations 2 and 14. For the original form of the equations, the 
evolution equation is 

V f ^ + u VB + BV ' u - B • Vu - uV • BJ = 0 

J^(V-B) = 0. (1C) 

l From the partial differential equation point of view, this might well seem the correct result; 
Y • B = 0 is an initial condition, and this equation guarantees that Y • B = 0 is maintained 
throughout the evolution. For the modified form of the equations, the evolution equation 
for the magnetic field is 

V- 1^ + u- VB+BV-u-B-VuJ = 0 

2- (Y • B) + V • (uY • B) = 0. (17) 
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Thus the addition of the source term has modified the evolution equation for V • B so that 
the quantity V • B jp is treated as a passive scalar. This is clearly the more numerically 
stable of the two evolution equations; any local V • B that is created is converted away. 

The above derivation gives all the pieces for building an ideal MHD solver that works 
for two-dimensional problems, without having to resort to non-collocated variables or a 
projection algorithm. Specifically, a Roe-type approximate Riemann solver has been imple- 
mented, where the wave strengths and speeds are derived from the above left eigenvectors 
and eigenvalues. The eigenvectors are properly normalized to avoid difficulties associated 
with coinciding wave speeds [5]. The average state needed at cell interfaces is computed by a 
simple average of left and right states (although a Roe average does exist for the ideal MHD 
equations [4]). The source term, though small, is calculated in each cell, and added to the 
residual. The resulting code is first order in space and time. 

3 A Test of the Eight- Wave Riemann Solver 

Brio and Wu [2] developed a test problem for one-dimensional MHD solvers based on the 
shock-tube problem of Sod [6]. Two stationary plasmas are separated by a membrane which 
is removed at t — 0, allowing the plasmas to interact. The test problem used here for the 
two-dimensional MHD solver is a rotated version of the Brio-Wu problem. The left and 
right input states, and the orientation of propagation of disturbances to the grid, is shown 
in Figure 2. In the Brio-Wu problem (the top figure), the boundary conditions are that the 
problem is periodic along a line y = constant; in the current test problem (the bottom figure), 
the boundary conditions are that the problem is periodic along a line x + y = constant. 

Both the rotated and non-rotated problems were run on coarse (COO cells in x) and fine 
(1200 cells in x) grids. The time step was taken as At./ Ax = 0.2, which corresponds to a 
OFL number of approximately 0.8 on the non-rotated problem. The ratio of specific heats, 
7, was 2.0. The number of time steps taken on the coarse and fine grids were 100 and 200, 
respectively. The x — axis in the plotted results from the rotated problem was scaled by a 
factor of \/2, to account for the fact that the CFL number is lower for the rotated problem 
than for the non-rotated problem. 

Figures 3-7 show comparisons of the results on the fine grid of the non-rotated shock-tube 
problem with the (scaled) results of the rotated shock-tube problem for 

3. density (p); 

4. pressure (p); 

5. velocity component normal to the original discontinuity (u„); 
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Figure 2: A Test Problem for Two-Dimensional MHD 
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Figure 3: Density in the Rotated and Non-Rotated Shock Tubes 
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Figure 4: Pressure in the Rotated and Non-Rotated Shock Tubes 
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Figure 5: Normal Velocity in the Rotated and Non-Rotated Shock Tubes 
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Figure 8: Normal Magnetic Field in tile Rotated Shock Tube (Coarse and Fine) 

6. velocity component tangential to the original discontinuity (u t ); 

7. magnetic-field component tangential to the original discontinuity ( B t ). 

As can be seen, the agreement is quite good, with the results of the two cases nearly indis- 
tinguishable for all but the normal component of velocity. The errors in u n are balanced by 
errors in the magnetic-field component normal to the original discontinuity ( B n ). Figure 8 
shows B n for the rotated shock-tube problem on the coarse and fine grids. In the non-rotated 
problem, B n = 0.75 throughout the tube. As can be seen, there are errors on the order of a 
few percent in B n on the coarse grid, but the errors are reduced as the grid is refined. 

4 Concluding Remarks 

In some respects, this paper presents the development of only one-eighth of a Riemann solver. 
Seven of the eight waves of the Riemann solver are the same as those used in previous work on 
upwind methods for MHD. The deceptively simple eighth wave that arises from the analysis, 
however, is of a different character than the other seven — it arises only in multi-dimensional 
problems, and it is crucial for understanding and solving those problems. It plays the very 
important role of stabilizing the numerical method with respect to the small amounts of 
V • B generated in solving the discrete MHD equations. 

Given the meteoric rise of Riemann solvers in the computation of compressible gas dy- 
namics, it is not very risky to predict that schemes based on Riemann solvers will play an 
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increasingly important role in the computation of compressible conducting flows. The ability 
of Riemann solvers to capture strong discontinuities robustly and with minimal dissipation, 
the framework that Riemann solvers provide for implementing stable boundary procedures, 
and the aesthetically appealing physical basis of Riemann solvers are all strong arguments 
for their use. The aim of this paper is to remove what is hopefully one of the last major 
obstacles to the use of Riemann solvers in large-scale codes for computing multi-dimensional 
conducting flows. 
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